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We study theoretically the effects of short-range electron-electron interactions on the electronic 
structure of graphene, in the presence of single substitutional impurities. Our computational ap¬ 
proach is based on the tt orbital tight-binding approximation for graphene, with the electron-electron 
interactions treated self-consistently at the level of the mean-field Hubbard model. We compare 
explicitly non-interacting and interacting cases with varying interaction strength and impurity po¬ 
tential strength. We focus in particular on the interaction-induced modifications in the local density 
of states around the impurity, which is a quantity that can be directly probed by scanning tunneling 
spectroscopy of doped graphene. We find that the resonant character of the impurity states near 
the Fermi level is enhanced by the interactions. Furthermore, the size of the energy gap, which 
opens up at high-symmetry points of the Brillouin zone of the supercell upon doping, is significantly 
affected by the interactions. The details of this effect depend subtly on the supercell geometry. We 
use a perturbative model to explain these features and find quantitative agreement with numerical 
results. 


PACS numbers: 73.22.Pr,71.55.-i,31.15.aq,71.10.Fd 

I. INTRODUCTION 

Graphene - a two-dimensional allotrope of carbon, has 
attracted considerable attention in recent years, largely 
due to its remarkable electronic properties stemming 
from the massless-Dirac-fermion nature of its low-energy 
quasiparticle states.® Most of the electronic properties of 
graphene that have been studied experimentally can be 
well described by non-interacting single-particle theory. 
However, electron-electron interactions in graphene are 
expected to be strong. In undoped clean graphene, the 
density of states at the Fermi level vanishes and therefore 
the Coulomb potential is not screened.®^ Recent experi¬ 
ments have shown that unscreened Coulomb interactions 
lead to reshaping of the ideal conical energy dispersion 
expected in graphene.® More precisely, the Fermi velocity 
near the Dirac point acquires a logarithmic correction as 
a result of interactions. 

From a theoretical viewpoint, it can be shown that 
this logarithmic enhancement arises from the non-local 
exchange interaction, already at the level of the first- 
order Hartree-Fock perturbation theory.® Hence, it is the 
long-range nature of the electron-electron interactions in 
graphene that is responsible for the logarithmic correc¬ 
tion, the most striking interaction effect observed so far 
in this material in the absence of external magnetic fields. 
As a result, theoretical work has mostly focused on inves¬ 
tigations of long-range interactions in graphene, using a 
variety of techniques ranging fro m mean fielcP to renor¬ 
malization group approaches.®^ 

It should be noted, however, that there are several im¬ 
portant conditions that need to be satisfied in order to 
observe significant long-range interaction effects exper¬ 
imentally. It is necessary to be able to probe a wide 
range of carrier concentrations and to tune the Fermi 
level sufficiently close the Dirac point, where the renor¬ 
malization of the Fermi velocity is expected to be dra¬ 
matic due to the vanishing density of states. Moreover, 


spurious screening effects, e.g. dielectric screening from 
the substrate, should be avoided. This makes undoped 
high-quality suspended graphene an ideal platform for 
studying the effects of long-range electron-electron inter¬ 
actions.® On the contrary, in the case of graphene on a 
substrate or in the presence of disorder and impurities 
these effects are less relevant. 

In particular, doping introduces a finite density of 
states at the Dirac point of graphene and the long-range 
part of the Coulomb potential is screened. In this case, 
short-range interactions become crucial. If these interac¬ 
tions are fairly strong, they can lead to interesting effects 
on the electronic structure, especially on the impurity 
states in the vicinity of the Fermi level. In fact, estimates 
of the o n-sit e Hubbard U parameter in carbon-based 
molecule^^Iini suggest that the short-range Coulomb in¬ 
teractions among 7r-electron in graphene can be indeed 
quite large, i.e. of the order of 10 eV. A similar value is 
obtained by accurate ab initio calculations.^^ 

In this paper, we study the effects of short-range 
electron-electron interactions on the electronic structure 
of doped graphene. We should note that the impor¬ 
tance of impurity effects in graphene has been addressed 
in many theoretical studies.C^Hill A number of interest¬ 
ing features have been revea led, including the opening 
of the gap upon dopinJ^^Ull ^nd the appearance of im¬ 
purity (acc eptor or donor) states in the vicinity of the 
Fermi level There is also a great interest in ad¬ 

dressing these properties experimentallj^MSI] since dop¬ 
ing graphene with impurities is one way to further explore 
and tune its electronic, magnetic and transport prop¬ 
erties. However, the interplay of short-range electron- 
electron interactions and impurity potentials in graphene 
has not yet been fully explored. 

We use a single-band (tt orbital) tight-binding (TB) 
model to describe the electronic structure of graphene. 
A supercell method is employed to study the effects of 
finite doping. A substitutional impurity is introduced 
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in the TB Hamiltonian as a local modification of the 
on-site potential at the impurity site. Here we focus on 
attractive impurity potentials, mimicking nitrogen impu¬ 
rity atoms which are typical dopants in graphene.^ The 
short-range interactions are described by means of the 
Hubbard model in the mean-field approximation, which 
is the simplest way of treating the many-body interacting 
problem. Interaction terms are introduced at each site in 
the TB Hamiltonian. We use a numerical self-consistent 
scheme to account for a redistribution of electronic charge 
around the impurity caused by interactions. As the out¬ 
put of our numerical calculations, we obtain the band 
structure and the local density of states (LDOS) around 
the impurity site. Furthermore, we calculate scanning 
tunneling microscopy (STM) images by integrating the 
LDOS over a small energy window above the Fermi level. 

By using this approach, we show that short-range in¬ 
teractions introduce several remarkable features in the 
electronic structure of doped graphene. Importantly, 
they enhance the resonant character of states localized 
in real space around the impurity, which are induced in 
the vicinity of the Dirac point. The complex interplay 
between short-range interactions and impurity potential 
is also responsible for non-trivial gaps at high-symmetry 
crossing points in the band structure of graphene, in par¬ 
ticular at the Dirac point. 

The paper is organized as follows. In Sec. |ll] we in¬ 
troduce our TB model and describe how the impurity 
potential and short-range interactions are incorporated 
in the Hamiltonian. We also provide some details of the 
self-consistent supercell calculations. Our findings are 
described in Sec. Ell In particular, in Sec. |HI A| we fo¬ 
cus on the effects of interactions on the band structure 
of graphene and on some issues related to the super¬ 
cell geometry. In Sec. |HIB| we discuss the changes in 
the resonant character of the LDOS around the impurity 
for varying impurity potential strength and interaction 
strength. The comparison between the simulated STM 
topographies for non-interacting and interacting cases is 
provided. Finally, we draw some conclusions. 


II. METHODOLOGY 

The second-quantized Hamiltonian for interacting elec¬ 
trons on a honeycomb lattice in the presence of impurity 
can be written as 

i(7 ^ 

( 1 ) 

Here cj^ and Cia are the creation and annihilation op¬ 
erators for electron on site i and with spin cr; Si and 
tij are on-site energies and hopping parameters, respec¬ 
tively. Only hopping between nearest neighbors on the 
honeycomb lattice is included. We assume that the TB 
parameters are uniform, except for the on-site energy at 
the impurity site, and we use the values obtained by fit¬ 


ting the TB band structures to density functional theory 
calculations, namely Si = 0 and tij = —2.97 

The third term in Eq. 0 represents the local impurity 
potential, with 17;^ being the impurity potential strength 
(Uim < 0 for attractive impurity). In our calculations we 
use Uim = —10 eV and Him = —20 eV in order to obtain 
visible trends for the impurity states in the vicinity of 
the Fermi level. 

The last term describes the on-site interaction between 
two electrons with opposite spins on site i (including the 
impurity site), with U {U > 0) being the Hubbard U pa¬ 
rameter, which expresses the strength of the intra-atomic 
Coulomb repulsion. Here riia- is the number operator, 
defined as riia = cl^Ci^. We consider H = 0, or non¬ 
interacting case, and U = 9.3 eV, which is the value ob¬ 
tained for graphene using the constrained Random Phase 
Approximation method.^ In order to extract the trends 
in the electronic structure with increasing the interaction 
strength we also use a larger value of H = 20 eV. 

In the mean-field approximation, the two-body inter¬ 
action term in Eq. Q becomes 

= U^(^{nii)cljC^j + (n,^)4^C4), (2) 

i i 

where (riia-) is the average electron occupation number, 
or density, for spin-up (u =t) and spin-down (cr =),) 
electrons. Here we consider a non spin-polarized case so 
that (riij) = (nil). 

In pristine graphene with the Fermi level exactly at the 
Dirac point, the average electron occupation number is a 
constant equal to 1/2. Adding a mean-field field on-site 
potential does not break the translational invariance of 
the crystal and the average occupation number remains 
constant. In fact, in orthogonal basis such a potential 
merely introduces a rigid shift of the energy bands (note 
that in non-orthogonal basis the interplay between the 
overlap integrals and the on-site interactions can lead to 
renormalization of the Fermi velocitj®l). 

However, the presence of both mean-field on-site in¬ 
teractions and impurity potential can lead to non-trivial 
effects in the electronic structure. In this case, the po¬ 
tential at each site depends on the average occupation 
number (riia), which is not necessarily the same on all 
sites. As a result, when a carbon atom is replaced by 
an impurity, there will be a redistribution of electronic 
charge in the system. In order to capture this effect, 
we need to perform self-consistent calculations for the 
Hamiltonian in Eq. Q and Q. 

At each step of the self-consistent cycle, the average 
occupation number for site i is calculated as 

- occ 

(3) 

k 

where N is the number of fc-points in the Brillouin zone 
and bi^ are the coefficients in the expansions of the wave- 
functions of the Hamiltonian in terms of the localized 
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atomic orbitals jicr). These are obtained by diagonaliza- 
tion of the Hamiltonian at each fc-point. The sum runs 
over all occupied states up to the Fermi level. Note that 
all calculations are done at half-filling. 

As initial values we use the occupation numbers calcu¬ 
lated for a non-interacting problem, i.e. for a supercell of 
graphene with impurity (C/im 7 ^ 0) and with U = 0. The 
criterion of self-consistency is 




(^zcr) 


s — 


1 


< V, 


( 4 ) 


where s is the index of the self-consistent cycle and 77 is a 
small parameter (we choose f] = 10“^). We use a linear 
mixing scheme, in which the input density (nicr)f„ at 
step s -|-1 is calculated as a linear combination of outputs 
(nicr)out ^'^‘4 (^io-)out^ from two previous steps 

= (1 - A) (ri*<T)out^ + A , (5) 

where A is the mixing coefficient; we use A = 0.25, 
which allows us to achieve self-consistency in less than 
100 steps. 

In order to model the effect of finite impurity con¬ 
centration, we construct a p x p supercell by replicat¬ 
ing a graphene unit cell p time along each of the two- 
dimensional lattice vectors [see Fig. [^a)]. The impurity 
atom substitutes a carbon atom in the supercell. In this 
work, we use two different supercells with p = 6 and 
p = 7. Atomic concentration of the dopants depends on 
the size of the supercell so the concentration is slightly 
different for the two choices, namely 1.0% for a 7 x 7 and 
1.4% for a 6 X 6 supercell. It is known that for p = 3q, 
where q is an integer, the Dirac points of graphene, K and 
K\ are mappe d onto t he F point of the Brillouin zone 
of the supercell.lA^HlIIIl 3,3 illustrated in Fig. S'*) . This 
does not happen if p is not divisible by 3. Therefore, the 
6 x 6 supercell is special. As we explain in Sec. 
the effects of impurity potential and interactions in this 
case are rather non-trivial. This is the main reason for 
considering two different supercell sizes. 



FIG. 1. (Color online) (a) 6 x 6 supercell of graphene with 
a substitutional impurity. A magenta sphere represents the 
impurity atom. Dashed lines mark the unit cell of pristine 
graphene and arrows show the primitive lattice vectors. A 
and B denote carbon atoms in the two equivalent sublattices, 
(b) Brillouin zone folding in graphene. Shaded area (1) repre¬ 
sents the first Brillouin zone of a p x p supercell of graphene. 
Numbered curves correspond to the first Brillouin zone of the 
unit cell for different p: p = 3g — 1 (2), p = Sq (3) and 
p = 3g -I- 1. For p = 3q the Dirac points of graphene {K, 
K') are mapped onto F point of the folded Brillouin zone. 
In other cases, K and K' are mapped onto K and K' of the 
folded Brillouin zone. 


III. RESULTS 
A. Bandstructure 

It is known that the finite amount of dopin g opens 
up an energy gap at the Dirac point of graphene! ^^ l 45 | i 7 | 
Here we address the question of how the details of the 
bandstructure near the gap are affected by interactions. 

We start with a special supercell geometry pxp, with p 
divisible by 3 (p = 6 in our calculations). Figure [^a)-(b) 
shows the bandstructure of the 6 x 6 supercell with impu¬ 
rity potential Him = —10 eV and Him = —20 eV, respec¬ 
tively, for three values of the interaction strength, H = 0, 
H = 9.3 eV and H = 20 eV. Note that in the bandstruc¬ 
ture calculations, different impurity potential and inter¬ 
action strength introduce a shift of the energy bands with 


respect to a reference case, i.e. non-interacting pristine 
graphene (Him = 0 and H = 0). In order to examine the 
features around the gap for different choices of parame¬ 
ters, we align the position of the doubly degenerate state 
(see the discussion below) in all curves in Fig.[^a)-(b) to 
the value found for H = 0 for a given impurity potential 
strength. 

As we mentioned before, in th e case of p = 6 both K 
and K' are mapped onto F point ,l42lM25| producing four 
degenerate states at F in the absence of impurities and 
interactions. When one carbon atom in the supercell is 
substituted by an impurity atom, a gap opens up between 
two states at F, however the other two states remain 
degenerate. More precisely, for H = 0 three of the four 
states at F are degenerate while one of the states moves 
away from the Dirac point. This situation is referred to 
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(a) 6X6 a^-lOeV (b) 6X6q„=-20eV 



FIG. 2. (Color online) Bandstructure of doped graphene 
along the high-symmetry lines of the Brillouin zone for 6x6 
(a,b) and 7x7 (c,d) supercell, for varying impurity poten¬ 
tial strength f/im and interaction strength U. Left panels are 
for Uim = —10 eV, right panels for t/im = —20 eV. In each 
panel three cases are shown: [7 = 0 (black), U = 9.3 eV (red) 
and [/ = 20 eV (blue). Horizontal line in (a) and (b) is the 
position of the doubly degenerate state at F, adjusted to the 
value found for [7 = 0 (see text for details). Horizontal line 
in (c) and (d) marks the conduction band maximum, which 
has been aligned with the value found for U — 0. 


as the pseudoga'^^ since there is still a pair of linearly 
dispersed states crossing at the neutrality point. Hence, 
effectively there is no band gap for this special supercell 
size. 

One can clearly see from Fig. Sa) that the size of 
the pseudogap decreases with increasing the interaction 
strength. On-site interactions cause a redistribution of 
charge around the impurity. We find that in the case of 
the attractive impurity potential Uim = —10 eV, the total 
average occupation at the impurity site decreases from 
(no) = 0.89 in the non-interacting case to (no) = 0.82 
for U = 9.3 eV. The on-site Coulomb repulsion pre¬ 
vents extra charge from accumulating on the impurity 
and therefore the strength of the impurity potential is 
effectively reduced. For the impurity potential which is 
twice stronger, the pseudogap for the same values of U 
is noticeably larger [see Fig.j^b)]. 

In addition to a large pseudogap, for 17 yl 0 there is also 
a smaller pseudogap which opens up at F [see the insets 
in Fig. [^a)-(b)]. There are now only two generate states 
at F, while the other two states shift, respectively, above 
and below the crossing point. Interestingly, the effect of 
interactions on the small pseudogap is opposite to that 
on the large one, e.g. its value increases with increasing 
the interaction strength. A perturbative model described 
in Appendix [A| suggests that the smaller gap results from 
the contribution of the states localized on sublattice B, if 
we assume that the impurity is substituted in sublattice 


A. Within our model, this effect is solely due to interac¬ 
tions (the contribution of sublattice B to the gap is zero 
in the absence of interactionJi^. We find the following 
values for the two gaps at F from analytical calculations 
(see Appendix [A] for details). In the non-interacting case 
and Uim = —10 eV, the large pseudogap is -0.56 eV. For 
U = 9.3 eV, the large pseudogap decreases to 0.39 eV. At 
the same time, a small pseudogap of 0.05 eV opens up. 
For U = 20 eV, the large and small pseudogaps become 
0.28 eV and 0.06 eV, respectively. These values are all in 
good agreement with the pseudogaps found in Fig. [^a). 
Note that for the larger impurity potential, the trends 
with increasing U are the same, however for a given U 
both gaps are larger than in Um = —10 eV case. Ana¬ 
lytical calculations using the perturbative model in this 
case also agree with numerical results. 

Similar features are found for a regular 7x7 supercell 
[see Fig. [^c) and (d)]. Note that the conduction band 
minima at F have been aligned with the reference U = 0 
case. The main difference from the 6x6 supercell is that 
in this case there is a real band gap at K {K'). Our 
perturbative analysis shows that there is no contribution 
from sublattice B to the gap at the Dirac point, if we 
assumed that the impurity is substituted in sublattice 
A. As in the case of the 6x6 supercell, the size of the 
gap decreases with increasing the interaction strength. 
Analytically, for Um = —10 eV we hnd a gap of 0.20 
eV for U = 0, 0.14 eV for U = 9.3 eV and 0.10 eV for 
77 = 20 eV. These values are in good agreement with 
numerical calculations. Somewhat smaller values of the 
gaps compared to a 6 x 6 supercell for the same Uim and U 
are expected since the atomic concentration of impurities 
is smaller. 

Bandstructure calculations presented in this section 
lead to a conclusion that short-range interactions ef¬ 
fectively reduce the strength of the impurity potential, 
which results in a decrease of the large gap (pseudogap) 
at the Dirac point. In oder to see how the character, 
e.g. the energy and the spatial extent, of the electronic 
states around the Dirac point is affected by interactions, 
we need to look at the LDOS around the impurity. 


B. Local density of states 

Calculations of LDOS at the impurity site reveal sev¬ 
eral important features. A substitutional impurity in¬ 
troduces electronic states at energies comparable to the 
impurity potential (jt/imj 10 eV), i.e. far away from the 
Fermi level. However, there are also states ap pearing in 
the vicinity (within ~ 1 eV) of the Fermi level.O^ ^ * ^^ * ^^ ! 
These states are the most relevant for the low-energy 
electronic properties of graphene and will be examined 
in detail. 

Figure [^shows the double- or multi-peak impurity res¬ 
onances close to the Fermi level for the 6x6 and 7x7 
supercells, respectively, for different impurity potential 
and interaction strengths. The multi-peak structure of 
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FIG. 3. (Color online) LDOS of doped graphene at the im¬ 
purity site for 6x6 (a,b) and 7x7 (c,d) supercell, for vary¬ 
ing impurity potential strength t/im and interaction strength 
U. Left panels are for Uim = —10 eV, right panels for 
Uiui = —20 eV. In each panel three cases are shown: [7 = 0 
(black), U = 9.3 eV (red) and t/ = 20 eV (blue). Vertical lines 
mark the position of the Fermi level (see text for details). 


the impurity resonances most likely originate from the 
long-range interaction, or interference, between the im¬ 
purity potentials, caused by the periodicity of the super¬ 
cell geometry.li^ 

With increasing the impurity potential strength, the 
resonant peaks move closer to low energies. This is very 
similar to the case of strong potential impurities on the 
surface of a topological insulator .1^ As shown in Fig.j^a) 
and (b) with [/ = 0, the double-peak resonance in the 
6x6 supercell approaches the Fermi level as Uim in¬ 
creases. At the same time its amplitude decreases. This 
finding is in agreement with semi-analytical calculations 
in Refs. miMl andUZl The effect of the impurity poten¬ 
tial is more striking in the case of a regular 7x7 supercell 
[Fig. I^c) and (d)]. In this case, in addition to shifting 
the peaks to lower energies, a stronger impurity poten¬ 
tial Uim = —20 eV also makes the peaks narrower, thus 
enhancing their resonant character [see in particular the 
first peak in Fig. [^c) and (d) with [7 = 0]. 

In the interacting case, the amplitude of the resonances 
increases with U. This can be seen in all panels in Fig. [^ 
with U = 9.3 eV and [7 = 20 eV, with the exception of 
the 6x6 supercell with Um = —10 eV and [7 = 20 eV, 
where the amplitude of the peak decreases slightly. In 
the case of Um = —10 eV, for both supercells the impu¬ 
rity resonances move further away from the Fermi level 
with increasing U. This is perfectly consistent with our 
observations for the non-interacting case with decreasing 
impurity potential. However, in the case of a very large 
impurity potential Um = —20 eV, the trend in the posi¬ 
tion of the resonances is less obvious. The peaks either 
do not move appreciably as in the case of [7 = 20 eV 
or even seem to move slightly towards the low-energy re¬ 


gion for [7 = 9.3 eV. Below we elaborate more on these 
findings. 

Increasing U reduces the overall strength of the im¬ 
purity potential, which is confirmed by decrease of the 
energy gap at the Dirac point due to the presence of 
impurities (Sec. Ill A). However, short-range electron- 
electron interactions controlled by U do not only change 
the potential directly at the impurity site but also affect 
the on-site potential and the charge density around the 
impurity (primarily nearest and next-nearest neighbors 
of the impurity atom). Hence, both the amplitude and 
the spatial extent of the impurity potential is modified 
by interactions. Let us assume that an attractive impu¬ 
rity can be described by a delta-function potential well. 
When interactions are included, the shape of the impurity 
potential is smoothed out (it acquires, say, a Gaussian 
shape). Therefore, although the strength of the poten¬ 
tial is reduced by a certain amount with increasing [7, 
the potential can become more long-ranged (in a certain 
parameter space). This, in turn, will increase the overlap 
of the potentials from neighboring cells and enhance the 
inter-supercell interaction. 


This seems to be the situation for Um = — 20 eV and 
[7 = 9.3 eV, for both choices of the supercell. In this case, 
the impurity potential decreases slightly due to interac¬ 
tions ([/ < Um), leading to a small decrease of the gap 
at the Dirac point compared to [7 = 0 case [Fig.j^b) and 
(d)]. At the same time, we find a significant difference be¬ 
tween the average occupation numbers of the nearest and 
next-nearest neighbors for [7 = 0 and U = 9.3 eV (this 
can be also partly seen in the STM images in Fig. [^a)- 
(b) and (d)-(e), for states in a small energy window above 
the Fermi level, where the neighbors of the impurity site 
appear brighter in the U = 9.3 eV case compared to 
[7 = 0). As a result, the amplitude of impurity res¬ 
onances in LDOS increases but their position shifts to 
lower energies. 

In contrast to this, for Um = — 20 eV and [7 = 20 eV, 
the average occupation numbers of the nearest and next- 
nearest neighbors of the impurity site do not change ap¬ 
preciably compared to [7 = 0 case. The strength of 
the impurity potential for this value of U is significantly 
reduced, leading to a large decrease of the energy gap 
[Fig. [^b) and (d)]. As a result, the amplitude of impu¬ 
rity resonances increases further, however their position 
remain close to the [7 = 0 value. These features strongly 
suggest that the position of the impurity resonances is 
sensitive to the spatial extent of the impurity potential, 
namely the resonances move closer to zero energy when 
the potential becomes more long-ranged. 

To further clarify the changes in the intensity and 
the spatial character of the low-energy impurity peaks, 
brought about by interactions, we present the simulated 
STM topographies in Fig. j^and Fig. [^ for Um = —10 eV 
and Um = —20 eV, respectively. For this we plot LDOS 
for each atom in the supercell,^ integrated over the 
energy window AE above the Fermi level (we choose 
AE = 0.25 eV). This gives an estimate of the tunnel- 
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FIG. 4. (Color online) Simulated STM topographies (LDOS 
for all atoms in the supercell) for 6x6 (left panels) and 7x7 
(right panels) supercell, for a fixed impurity potential strength 
Uira = —10 eV and varying interaction strength: [7 = 0 (a,d), 
U = 9.3 eV (b,e) and 1/ = 20 eV (c,f). Arrows mark the 
position of the impurity atom. Note the logarithmic color 
scale. 


ing current as electrons tunnel out of the occupied states 
of the STM tip into the unoccupied states of graphene. 
The increase of the electronic density of states in this 
energy window gives rise to a bright triangular feature 
around the impurity. Note that for the 6x6 supercell 
the impurity is located in sublattice A, while for the 7x7 
supercell it is located in sublattice B. Therefore the bright 
triangular features in the two supercells appear rotated 
by 180° with respect to each other. 

The difference between non-interacting and interact¬ 
ing cases is clearly visible in the STM images. In all 
cases considered, the impurity site becomes progressively 
brighter compared to its neighbors with increasing U. 
This means that the electronic states in the small energy 
window above the Fermi level become more localized on 
the impurity site as a result of interactions. This is most 
evident in the case of the 7x7 supercell [Fig. [^d)-(f)]. 
For C/im = —10 eV, the overall intensity of the images 
decreases with U. For a stronger C/im = —20 eV impu¬ 
rity potential, the trend is similar with the exception of 
the intermediate interaction strength U = 9.3 eV. In this 
case, the contribution of the nearest and next-nearest 



FIG. 5. (Color online) Same as Fig. but with impurity 
potential strength Uim = —20 eV. 


neighbors is even stronger than in the U = 0 case. This 
correlates with the corresponding features in the LDOS 
discussed above. 


IV. CONCLUSIONS 

We have presented a theoretical study of the effects 
of electron-electron interactions on the electronic states 
of graphene in the presence of substitutional impurities. 
Using a self-consistent TB model with on-site interac¬ 
tions treated at the mean-field level, we have shown that 
the size of the gap, which opens up at the Dirac point 
in graphene upon doping, and the character of the low- 
energy electronic states are modified by interactions. The 
mechanism for these effects is provided by the interplay 
between the impurity potential and the on-site repul¬ 
sion, which leads to significant re-arrangement of the 
electronic charge around the impurity compared to the 
non-interacting case. 

In particular, we found that the size of the gap de¬ 
creases with increasing the interaction strength. Intu¬ 
itively, this can be understood as follows. In the case of 
an attractive impurity potential, which mimics nitrogen 
dopants in graphene, adding the on-site Coulomb repul¬ 
sion effectively reduces the strength of the potential, i.e. 
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the depth of the potential well decreases. This is due 
to the fact the on-site repulsion prevents extra electronic 
charge from accumulating on the impurity site. 


For a special supercell size pxp, where p is divisible by 
3, both K and K' are mapped onto F point of the folded 
Brillouin zone. Therefore, in the case of undoped non¬ 
interacting graphene, there are four degenerate states at 
the neutrality point. It is known that when the impu¬ 
rity potential is included, a gap (pseudogap) opens up 
between two of these states while the other pair remains 
degenerate.^ We have shown that the size of this pseudo¬ 
gap is reduced by interactions. Interestingly, in addition 
to this, a small gap opens up between the second pair 
of states at the F point, which are otherwise degenerate 
in the absence of interactions. We explain these features 
both qualitatively and quantitatively, using a perturba¬ 
tive model based on the generalization of the approach 
developed by Lambin et to the interacting case. 

Furthermore, we have studied the behavior of the 
impurity-induced electronic states with and without in¬ 
teractions. There are two groups of states which can be 
detected in the density of states when a carbon atom 
in the supercell is replaced by an impurity atom. First, 
there are states which emerge far away from the Fermi 
energy, with their energies of the same order of mag¬ 
nitude as the impurity potential ( « 10 eV in our cal¬ 
culations). Second, there are states appearing close to 
the Fermi energy and are therefore of particular inter¬ 
est. Although the way the electron-electron interactions 
affect the LDOS at low energies in general depends on 
the impurity concentration, we found clear trends in the 
behavior of the impurity-resonances as both parameters, 
i.e. the interaction strength and the impurity potential 
strength, are modified. 


Regardless of the interactions, the impurity levels move 
closer to the low-energy region (e.g. to the original Dirac 
point) with increasing the impurity potential strength. 
This find ing is c onsistent with previous calculations for 
graphene.lI^EnilZl Similar result was found in another class 
of Dirac materials, namely in three-dimensional topolog¬ 
ical insulators in t he pre sence of strong potential impuri¬ 
ties on the surface.l^^*^ However, our self-consistent cal¬ 
culations for graphene in the presence of both impurities 
and interactions reveal novel features. For a fixed im¬ 
purity potential, sort-range interactions tend to enhance 
the amplitude of the impurity-resonances in the vicin¬ 
ity of the Fermi level. The position of the resonances is 
also affected by the spatial extent of the effective impu¬ 
rity potential, which is modified by interactions. As the 
interaction strength increases, the states become more 
localized on the impurity atom. The differences in the 
spatial distribution of the low-energy impurity states in 
the non-interacting and interacting cases are clearly de¬ 
tectable in the simulated STM topographies. 


Appendix A: Band gap in doped graphene supercell 
in the presence of interactions 

We use a simple perturbative model, based on the 
model proposed in Ref. IT3] for non-interacting graphene, 
to explain how the impurity potential and electron- 
electron interactions affect the band structure of doped 
graphene. 

We first review this model for the non-interacting case 
with impurity.^ The Hamiltonian for such a system can 
be written in Dirac’s notation as 

H = Ho + H,, (Al) 

where Hq is the Hamiltonian of non-interacting pristine 
graphene 


Ho = Y,HE^{u\+Y, |w) tuv (uj, (-^2) 

u u,v 

and Hi is the perturbation introduced by the periodic 
arrangement of impurities 

Hi = Y, \u) Uim {u\ . (A3) 

ItGl 


Here |u) is the atomic orbital associated with site u 
{{u\u') = Suu'), Eu and tuv are the on-site energies and 
the nearest-neighbors hopping integrals, respectively; 
t/im is the impurity potential and the sum over it G 1 
refers to all impurity atoms belonging to sublattice 1, 
which substitute one carbon atom in each px p supercell 


In pristine graphene [Eq. ( |A2[ )], there are four states 
with zero energy, two for the two sublattices and two 
for the non-equivalent points K and K' in the Brillouin 
zone (we omit the spin indices for simplicity). The cor¬ 
responding Bloch functions can be written as 


y e 


zK-i 


|u), (A4) 


where K is a vector in the reciprocal space corresponding 
to either A or AT', u is the position of site u in real space, 
and AfA(B) is the number of atoms in sublattice A(B). 

The task is to calculate the first order corrections to 
the energy states at the Dirac point due to the impurity 
potential, by using degenerate state perturbation theory. 
Let us assume that the impurity is substituted in sub¬ 
lattice A. Then the states and have zero 

amplitudes on atoms in sublattice A and these states are 
eigenstates of zero energy. Therefore, we only need to 
consider the subspace of degenerate eigenstate s fo rmed 
by t he s tates \K^) and \K'^). We use Eq. (A3) and 
Eq. (A41 to calculate the following matrix elements 


Vii = {K^\Hi\K^) = ^ = (A5) 

^11/ _/Va 

V22 = {K'^\Hi\K'^)=Vii = ^, (A6) 

Nil/ pZ 











Vi2 = (K^j Hi \K'^) = — 

uG 1 

= —^^K'-K,G, (A7) 

pz 

V21 = (K'^l Hi \K^) = V12 = %( 5 k'-k,g, (A 8 ) 

where Na = since for pxp supercell, we have N = 2p^ 
atoms and Na = Nq = p^ atoms in each sublattice. 
^ 12 (^ 21 ) in Eq. is not zero only when the vector 

K' — K is equal to the reciprocal lattice vector G. This 
condition is satisfied only if p is divisible by 3, i.e. if p = 
Sg, where p, q are integers. This is essentially the result 
obtained in Ref.[T3]and it is confirmed by our calculations 
presented in Fig. for [7 = 0 case. 

The first corder corrections eH) to the energy states 
are then given by the eigenvalues of matrix E, with 
matrix elements Vij {i,j = 1,2) defined above. This 
gives E^^'> = 0 and = 2Uim/p'^ if p = 3q, and 
= Uimlp"^ if P 7^ 3g. Hence, in the case when p 
is divisible by 3, all four states are mapped onto T point. 
Three of these states, two of which are localized on sub¬ 
lattice B and one on sublattice A, have zero energy. The 
remaining state is shifted down in energy by 2Ui^/p^ 
(Uini < 0 for attractive impurity), producing a pseudo¬ 
gap of magnitude 2?7ini/p^. This is exactly the situation 
shown in Fig. [^a) and (b) for the 6x6 supercell. In the 
case when p is not divisible by 3, the degeneracy between 
K and K' is lifted and there are now two states at each 
of these points, separated by a gap of C/im/p^. This is the 
situation for the 7x7 supercell in Fig. [^c) and (d). 

Combining the results for p = 3q and p ^ 3q, to the 
lowest order in the perturbation theory the gap (pseudo¬ 
gap) induced by the impurity potential can be written 
as 


Fgap = —^ -I- ^dK'-K,G- (A9) 

p2 pZ 

For Uiyn = —10 eV, Eg^p = —0.56 eV if p = 6 and Fgap = 
—0.2 eV if p = 7. These values coincide with the gaps 
found numerically [see Fig. ia) and (c) with [/ = 0]. 

Below we extend this model to the interacting case. 
The Hamiltonian of the interacting system can be written 
as 


H = Ho + Hi+H 2 


(AlO) 


where Hq and Hi are given by Eq. (A2) and Eq. 
respectively, and H 2 is the perturbation introduced by 
short-range interactions which is given by 


H2 = Y,\'^)iU-M){u\. (All) 

U 

In accordance with Eq. ([^, {riuo) is the average elec¬ 
tron occupation number on site u corresponding to spin 
a. When the mean-field interaction term is included, the 
energy bands acquire a rigid shift. In order to compare 


the results to the non-interacting case, we need to sub¬ 
tract this shift. Let us assume that it is proportional to 
an average quantity (h„o-)) which is a constant. Then the 
energy that needs to be subtracted is U■ {fiua)- 

We now calculate the first-order corrections to the en¬ 
ergy states at the Dirac point due to interact ions u sing 
the same procedure. Since the sum in Eq. (All) can 
be decomposed into the sum over u S A and the sum 
over u S B, we can calculate corrections due to these two 
terms separately. For each of them we need to consider 
either the subspace formed by the states and K'^, or 
by the states and K'^. Let us assume for simplicity 
that all occupation numbers for atoms in sublattice A 
are approximately the same and equal to {fiua-), except 
for the impurity site. Then the corresponding matrix 
elements are give by 


Wi\ = {K^\H 2 \K^) = ^ E ^ 


liG A,f7 

pZ 

pZ 

= {K'^\H2\ K'^) = W^i , (A14) 

Wi^i = {K'^\H 2 \K^) = IFi4, (A15) 


where u' is the impurity site; 2 in the denominator stands 
for spin. As before, first order corrections to the energy 
states are given by the eigenvalues of matrix VF^, with 
matrix elements ITA [j^j = 1,2), and differ for p = 3q 
and p ^ 3q. These corrections give the following contri¬ 
bution to the energy gap at the Dirac point 


U 


Agap^^ = ^ {{riu'a) — {nua)) (1 + I^K'-K.g) ■ (A16) 


In a similar way, we calculate the first-order corrections 
in the subspace formed by states and K'^. The cor¬ 
responding matrix elements are given by 

Wfi = {K^\H2 \k^) = ^ E - («-)) 

= ^ ^ ^ ((^UO-) ~ {fiua)) ) (A17) 

uCnn of u' 

Wfa = {K^\H2 \K''^) 

= ^ E {{nua} — {nua)) Sk'-K,G, (A18) 

P ^ t , 
nGnn oi u 

^22 = {K'^\H2\ K'^ ) = IFn, (A19) 

W^i = {K''^\H2\K^) = IFi®2 , (A20) 

where we assumed that all atoms in sublattice B have 
approximately the same occupation (fiua), except for the 
nearest neighbors of the impurity atom. This is a rea¬ 
sonable assumption since these atoms are most strongly 
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TABLE I. Values of band shifts and average occupation num¬ 
bers used in Eq. (A22l for Ui-m = —10 eV and U = 9.3 eV. 


Band shift 

(^4icr) 

(fiu'a) (imp.) 

(nua) (nn of imp.) 

4.60 eV 

0.49 

0.82 

0.46 


affected by the impurity. Then the contribution to the 
energy gap, stemming from corrections to the states lo¬ 
calized on sublattice B, is given by 


^ E (1 + ^K'-K,g) . 

liGnn of u' 

(A21) 

Finally, combining the corrections due to the impurity 
potential [Eq. (A9)] and due to interactions [Eq. (A16) 
and Eq. (A21)], we obtain the expression for the energy 
gap at the Dirac point 


E- 


gap 



(1 + ^K'-K.a)! 

{p^ua )) (1 + i5k'-K,g)| , 
(A22) 


where the expression inside the first curly bracket is due 
to the states localized on sublattice A, while the second 
one is due the states localized on sublattice B. 


Let us now summarize what happens to the four zero 
energy states, when both the impurity potential and in¬ 
teractions are present. When p = 3g, all four stat es ar e 
mapped onto F point. As one can see from Eq. (A22|, 


one of the states, corresponding to sublattice A, remains 
at zero energy, while the other one shifts by = 


pT + pi ' {{nu'a/ 


We can estimate this 


t) {P’ucr^^ 

quantity by taking the values of the rigid band shift 
and the average occupation numbers for impurity and its 
nearest neighbors from our numerical calculations (see 
Table |lj numbers are similar for the two supercells). For 
C/im = —10 eV, U = 9.3 eV and p = 6 , this energy shift 
is negative and is equal to —0.39 eV. This corresponds 
to the large pseudogap (below the Dirac point) that we 
identified in Fig. [^a) for this choice of U. In a similar 
way, one of the states localized on sublattice B remains 
at zero energy, while the other one shifts up in energy 

by Agap ^ = 2^ X^uGnn of u' — {fiua)) = 0.054 eV. 

This small energy shift is identical to the small pseudogap 
(above the Dirac point) found in Fig. [^a). 

When p ^ 3, the degeneracy between K and K' 
is lifted. There are two states at each of these 
points, one localized on sublattice A and the other one 
on sublattice B. The energy gap between the states 
is given by = ^ + E . ((n„/<,) - {fiua)) + 

^ X^uGnn of «' — (h-ucr))- For Dim = —10 eV, U = 

9.3 eV and p = 7, E“p = 0.14 eV, which is in good 
agreement with Fig. [^c). 
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